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An ever broader range of physical platforms provides the possibility to study and engineer quan¬ 
tum dynamics under continuous measurements. In many experimental arrangements the system of 
interest is monitored by means of an ancillary device, whose sole purpose is to transduce the signal 
from the system to the measurement apparatus. Here, we present a method of adiabatic elimination 
when the transducer consists of an arbitrary number of bosonic modes with Gaussian dynamics 
while the measured object can be any quantum system. Crucially, our approach can cope with the 
highly relevant case of finite temperature of the transducer, which is not easily achieved with other 
methods. We show that this approach provides a significant improvement in the readout of super¬ 
conducting qubits in circuit QED already for a few thermal excitations, and admits to adiabatically 
eliminate optomechanical transducers. 


I. INTRODUCTION 

Quantum limited, continuous measurements and 
measurement-based feedback [1, 2] represent important 
concepts for fundamental studies of open quantum sys¬ 
tems and the measurement process in quantum mechan¬ 
ics, and beyond that they are highly useful tools for ap¬ 
plications in Quantum Information Processing. Since the 
first demonstration of quantum limited, continuous mea¬ 
surements with single ions [3] and atoms [4] the concepts 
of quantum dynamics under continuous monitoring have 
gained great experimental relevance in recent years in the 
field of circuit QED [5] . Here, measurement and feedback 
have been used for preparation of qubit states [6-8] in¬ 
cluding preparation of entangled states [9, 10], or for ob¬ 
serving and stabilization of quantum trajectories [11-13]. 
Only very recently cavity optomechanical systems [14] 
entered the parameter regime of quantum limited, con¬ 
tinuous measurement [15] and feedback within the ther¬ 
mal decoherence time [16], where the tools of continuous 
measurement theory will unfold their full strength. 

In a typical measurement scenario, the system of inter¬ 
est (e.g., a superconducting qubit or a mechanical oscilla¬ 
tor) interacts with an ancillary system, such as a cavity 
mode, whose continuously emitted output field is then 
sent to a measurement device, cf. Fig. 1. In the most 
simple case the ancilla is a single cavity mode transduc¬ 
ing the signal photons emitted from e.g. a superconduct¬ 
ing qubit to the measurement apparatus. Beyond that, 
the ancillary transducer can consist of a much more so¬ 
phisticated subsystem, e.g. an arrangement for frequency 
conversion of photons from microwave to optics, which is 
subject to its own nontrivial dynamics, losses and sources 
of thermal noise. The precise dynamics of the auxiliary 
system is oftentimes irrelevant and we are interested only 
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Figure 1. (a) Schematic of the considered setup: A quantum 
system is monitored by coupling it to an ancillary system - 
the transducer - whose continuously emitted light field is de¬ 
tected. The most simple example of such a setup—a qubit 
coupled to a cavity mode with monitored output—is shown 
in (b). The transducer can also be a much more complex de¬ 
vice, e.g. an optomechanical converter of microwave to optical 
photons. 


in obtaining the equation of motion of the system alone. 
Obtaining the reduced dynamics of the system is cru¬ 
cial for two reasons: Firstly, auxiliary systems quickly 
make any numerical simulations intractable due to the 
increased Hilbert space dimension. This becomes espe¬ 
cially troublesome when the transducer is thermally ex¬ 
cited at finite temperatures. Secondly, in more complex 
setups, where the transducer is composed of several cou¬ 
pled subsystems, the structure of the ancilla easily ob¬ 
scures the effect of the continuous measurement on the 
main system, and makes it difficult to design feedback 
protocols. 

The dynamical degrees of freedom of the transducer 
can be adiabatically eliminated from the dynamics if their 
evolution (e.g. the cavity decay) is fast on the time 
scale of the interaction with the system. In the con¬ 
text of stochastic quantum dynamics under continuous 
measurements perturbative techniques used for adiabatic 
elimination of fast degrees of freedom received signifi¬ 
cant attention in the theoretical literature in this field 
[17-24]. However, these methods become cumbersome 
or intractable when the ancillary system becomes large 
(i.e. consists itself of several subsystems) and/or ther¬ 
mally excited at finite temperature. Both of these cases 
are highly relevant for current experiments: On the one 
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hand, thermal excitations typically cannot be neglected 
in the frequency domain of circuit quantum electrody¬ 
namics [5]. On the other hand, optomechanical systems 
composed of several coupled mechanical, optical and mi¬ 
crowave modes can be used as transducers in order to 
convert photons at vastly different length scales [25-30]. 

In this paper, we present a new method for adiabatic 
elimination in conditional dynamics which applies to 
transducers composed of an arbitrary number of bosonic 
modes whose dynamics is Gaussian [31, 32], that is, the 
dynamics is generated by a quadratic Hamiltonian, linear 
jump (decay) operators, and the transducer is subject to 
continuous homodyne detection. We use the fact that 
dynamics of such systems can be described using their 
first and second statistical moments—the mean values 
and the covariance matrix—to obtain a stochastic mas¬ 
ter equation for the system of interest only. The main 
advantage of this method is the possibility to eliminate 
a subsystem of arbitrary dimension coupled to thermal 
bathes which implies a broad range of applications for 
superconducting and optomechanical systems. 

We introduce the method in Sec. II, and illustrate the 
method with several examples in Sec. III. Not only do 
we show that our method of adiabatic elimination outper¬ 
forms more naive approaches when dispersively reading 
out a qubit using a cavity with just a few thermal excita¬ 
tions, we also study measurement-induced entanglement 
generation in such a setting, generalizing the results of 
Ref. [19]. Finally, in Sec. IV we conclude and suggest 
other applications of our method in quantum control sce¬ 
narios with optomechanical transducers [25-30]. 

II. EFFECTIVE DYNAMICS 
A. Main results 

The measurement scenario we have in mind is illus¬ 
trated in Fig. 1. The system of interest (e.g., a qubit) 
couples to a transducer (e.g., a cavity mode) whose out¬ 
put fields are continuously measured in a homodyne de¬ 
tection. The conditional dynamics of the overall system, 
including losses, noise, and the effect of continuous diffu¬ 
sive measurement, is described by the stochastic master 
equation 

dp = Cspdt + Crpdt + Cmtpdt ^ n[Xm]pdWm- (1) 

m 

Cs is the system Liouvillian that contains, in general, 
some coherent dynamics given by a system Hamiltonian 
Hs and some Lindblad operators decribing decoherence, 
but will be left unspecified for the moment. The Liouvil¬ 
lian for the transducer is 

CtP= -i[HT,p]+Y^V[ji]p, 

i 

where Hx is the Hamiltonian and the Lindblad terms 
T^ljilP = jiPj] - ^JiJiP + PJiJi) describe decoherence 


and measurement channels. We further assume that the 
transducer is Gaussian, i.e. it consists of N bosonic 
modes with Hamiltonian Ht and jump operators ji which 
are, respectively, quadratic and linear in canonical oper¬ 
ators. It is convenient to collect the canonical operators 
into the 2V-dimensional vector r = (gi,pi,..., qM^PN^ 
with commutation relations [ri,rj] = icrij and cr = 

^ being the symplectic matrix. The trans¬ 

ducer Hamiltonian and the jump operators can then be 
expressed as, respectively, 

Ht = 

with a real symmetric matrix R = G x and 
complex vectors Furthermore, we assume the 

interaction between system and transducer linear in the 
transducer operators 

^intP — p\j .^int — ST. 

where s is a 2iV-dimensional vector of Hermitian opera¬ 
tors acting on the system S. We use the small parameter 
e to remind us that the interaction is weak and can be 
treated perturbatively. Finally, for the transducer to be 
Gaussian, the measurement terms correspond to a homo¬ 
dyne detection, i.e., 

T~{-[Xm]p = (Xm ~ (Xm))p + P(Xl„ — (Xj^)), 

and the measurement operators are linear in the canoni¬ 
cal operators, 

Xjn — (Cm T r, (2) 

with CrmTUm G The measurments are independent, 

dWmdWn = Smndt and for each measurement opera¬ 
tor Am, there should be a corresponding Lindblad term 
D[Am] in the transducer Liouvillian Ct- The measure¬ 
ments give rise to classical measurement currents that 
take the form 

dim = (Am + A]j)j)dt -|- dWm- (3) 

To zeroth order in the coupling parameter e, the trans¬ 
ducer dynamics is Gaussian which means that it can be 
fully described using the first and second statistical mo¬ 
ments of the canonical operators, i.e., the mean values 
Xi{t) = {ri(t)) = tY{pT{t)ri}^ and the covariance matrix 
with an element I'ijit) = ([r^,rj] + (t)) — 2xi{t)xj{t) = 
tT{pT{t)[ri,rj\ + } - 2xi{t)xj{t)\ here px = trs{p} is the 
reduced density operator of the transducer (in zeroth or¬ 
der of e) and we use the superscript c to indicate that the 
moments are calculated with respect to the conditional 
state pt, obeying the stochastic master equation 

dpx = C-xprdt -I- ^ 'H[Am]pTdWm. (4) 
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For the mean values and the covariance matrix this im¬ 
plies the following equations of motion 

dx = AxAt -I- Cm - amm)dWm, (5a) 

m 

f^ = AT^ + + 2A^ - 

2(r Cm Cm 

where 

A = aR- (6a) 

i 

(6b) 

i 

see Appendix A for details. 

Averaging Eq. (4) over the measurement record, we 
recover the deterministic master equation for the uncon¬ 
ditional state 

p't ~ ^tPt 

(we use the superscript u to indicate that the state is un¬ 
conditional) and the corresponding equations of motion 
for the first and second moments 

= Ax'^, (7a) 

F" = AT'^ + -f 2N. (7b) 

Note that in both cases, conditional and unconditional 
dynamics, the covariance matrix obeys a determinis¬ 
tic equation of motion of Ricatti or Lyapunov type, cf. 
Eqs. (5b) and (7b) respectively, which can be solved efh- 
ciently. 

Our main goal is to derive a closed, effective equa¬ 
tion of motion for the conditional state of the system 
PS = trT{p} which is correct to leading order of e based 
on the assumption that transducer dynamics Ct is fast 
on the time scale of the system-transducer interaction 
Hint- Under this condition the state of the system will 
be given by p = ps 0 pr + 0{e) where pr is the steady 
state solution of Eq. (4). The strategy now is to deter¬ 
mine equations of motion for the order-e correction to 
this approximation, solve them formally, and substitute 
the solution into the equation of motion for ps- In this 
way we can arrive at the closed, effective equation of mo¬ 
tion for PS, which will be of second order in e in the 
deterministic and of first order in the stochastic part, as 
we will see. In the following we will summarize the fi¬ 
nal result of this adiabatic elimination procedure. The 
derivation is given in the next section. 

So far we have left the system dynamics Cs unspecified. 
Eor the adiabatic elimination to work we will have to 
make assumption regarding Cs relative to Hmt and Ct- 
We will consider two main regimes: 

(a) The system dynamics is trivial, £5 = 0 . This can 
be fulfilled exactly in an interaction picture when the 
system operators Sj in Hmt happen to be constants of 
motion, and covers in particular the important case of 


a quantum nondemolition measurement. Else, Cs = 0 
can be fulfilled approximately if the time scales of Cs 
are much slower than those of Hmt- Under this assump¬ 
tion the effective equation of motion for the state of the 
system ps is found to be 

1 i 

dps = [sfe,ps]]dt -k -Ay Vjfeh, [sfc,ps]-h]dt 

+ n[iAls]psdWm, ( 8 ) 

where we have for the measurement term 

Am = (r'" - ia)Q~'^Cm + A~^{T‘^Cm - crnim), (9a) 
Q = A-2{r^ Cm-cnrnm)Cm- (9b) 

We remind the reader that refer to the covariance 
matrix of the transducer attained as steady state solu¬ 
tions of Eqs. (5b) and (7b), respectively. The coefficients 
Cm and nim are given in (2), and the matrix A in (6). In 
the last equation and in following equations we use the 
Einstein summation convention. 

As expected, the deterministic part of the stochastic 
master equation (first line in Eq. (8)) depends only on 
the unconditional state of the transducer through its co- 
variance matrix F”. Note that the stochastic term does 
depend on the conditional state F'^. 

The effective equation of motion (8) is not manifestly 
in Lindblad form. In order to bring it into the Lindblad 
form we rewrite it as 

dps = -i[H,ps] + Pzj (^Sipssj - ^{sjSips + psSjSi)^ 

+ n[iA'^s]ps 

where 

H = (A“^(F'^ -k ia) — (F'^ — ia^)A~^) s, (10) 

P = -\ (A-1(F" - ia) + (F" + ta^)A-^) . (II) 

The individual jump operators and corresponding decay 
rates are given by eigenvectors Vi and eigenvalues Wi > 0 
of the matrix P, '^^Wi'D[vf s]ps- P is indeed a positive 
semidefinite matrix, as we show in App. B. 

Finally, the effective equation of motion has to be ap¬ 
pended with an equation relating the measured photocur¬ 
rent to the system observables Si after elimination of the 
transducer degrees of freedom (replacing Eq. (3)) 

dim = (*A^s - is^Al,}dt + dWm- 

(b) When the interaction and system Hamiltonians do 
not commute, moving to the interaction picture with re¬ 
spect to the system Liouvillian Cs results in a time de¬ 
pendent interaction. In the simplest and most common 
case the system operators oscillate at a particular fre¬ 
quency iw and the interaction Hamiltonian becomes 

= s'^{t)r, s{t) = s+e*“* -k 
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with time independent operators (s+)i = (s_)|. Since 
in this case the signal of the system (i.e. the pho¬ 
tons emitted by it) is now carried by sidebands it will 
be necessary to detune the local oscillators in the ho¬ 
modyne measurements, such that the canonical oper¬ 
ators in the measurement terms become time depen¬ 
dent, -|- a|e*'^”*‘)/-\/ 2 , pi = — 

aie“®^'"*)/-\/ 2 , with being the annihilation operator 
for mode i and the detuning of the local oscillator 
from the central frequency. The time dependence can be 
moved to the coefficients of the measurement operators 
in Eq. (1), Am = {cm + immYr{t) = {cmit) + r. 

Performing a coarse-graining in time (provided Am is 
faster than any other time scale of the transducer), the 
Riccati equation (5b) has to be replaced by 

pc = y4r= -h + 2N- (12a) 

m aG{c,s} 

Cm{t) = Cm cos(Amt) + sin(Amt), ( 12 b) 

mm{t) = TOm cos(Amt) + sin(Amt). ( 12 c) 

Eliminating the transducer, the system density operator 
then obeys the equation of motion 

dps = Cpsdt + 'H[Km]psdWm, (13a) 

d/m = (Am + + dWm, (13b) 

where the deterministic part is given by 

^Ps = -^{A + [s_,fc,ps]] -I- 

+* 0 ’jfc[S-|-,i, [S-.fc, Ps]-I-]) + 

+ 2 ~ [s+,k, Ps]] + 

+iajk[s-,z,[s+,k,ps]+]), (14) 

and the particular form of the measurement term de¬ 
pends on the choice of local oscillator detuning, for which 
one has to distinguish the two relevant cases Am = ±w. 


4” I Am — tU, 

(15a) 

4” I Am = W, 

(15b) 

(T'’ - io){Q + iAm)“'^c+ -f 


+ {A - iAm)~\T^c+ - crm+). 

(15c) 

(T'’ - ia){Q - iAm)“^Cm 4- 


+ {A + iAm)"^(r‘’Cm - CrrUm), 

(15d) 


Q = (r'=c^-crTO^)(c^)^, (15e) 

m aG{c,s} 

where we denote Cm(t) = 0 + 6 *^“* -I- €“ 6 “*^"** and 
are defined similarly. One should, once again, check that 
each measurement term has a corresponding Lindblad 
term in the unconditional part of the dynamics. 

In the rest of this section, we present detailed deriva¬ 
tions of the equations of motion Eqs. ( 8 ), (13). Reader 
interested in applications of these results may thus jump 


straight to Sec. Ill, where we illustrate the use of these 
equations on several examples concerning qubit readout 
in circuit QED. 


B. Adiabatic elimination with time-independent 
interaction 

We start the adiabatic elimination by simply tracing 
out the transducer dynamics from the stochastic master 
equation (I), leading to 

dps = trT(dp) = -ie[si, r]i]dt + 2cmiPidWm, (16) 
where we defined 

rji = irririp) - x^ps- (17) 

In view of pi = irTi^iP) — trT(»"iPT)trT{p} we can give 
a simple physical meaning to the quantities pp They 
measure the deviation of the exact state p from the ten¬ 
sor product state pr < 8 ) ps with respect to the first order 
moments of the transducer’s canonical variables r^. Ac¬ 
cordingly, for the tensor product state p = pr ® ps we 
have Pi = 0 , and, as we will see, the pi are of first order 
in e. Next, we derive equations governing the evolution 
of rji and pi, solve them formally to first order in e, and 
plug the solutions into Eq. (16). 

To obtain an equation for rji, we need to evaluate 

drji = tTrindp) 

= tTTinCTp)dt - ietTTiri[sjrj, p\)dt + 
+CmjtTT{ri[rj - Xj,p\ + )dWm + 
+immjtT:T{ri[rj,p])dWra 
2 

= AijTjj — —e(T^j + 2xiXj)\sj^ps\dt -\- 

+ ]^eaij[sj,ps\+dt -|- Cmj{Uij + 2xiPj)dWm - 
-aijrrimjPsdWm- (18) 

Here we used the fact that the first term on the right 
hand side of the above equation is completely analogous 
to terms appearing in the equation of motion for the 
mean values of the canonical operators Eq. (5a), and 
is therefore equal to AijPj. Eor the other deterministic 
part we further used 

= i([ri,rj]+ -h [n^rj]) 

~ 2 + icTij), (19) 

and p = PS ® Pt 0th order in e. Einally, we defined 
Uij = trx{[Ti — XijTj — Xj\+p). To solve Equ. (18), we 
also need an equation of motion for XiXjps that is valid 
to 0th order in e. Using the ltd product rule, d(Ay) = 
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{dX)Y + XdY + dXdY, we have 

dxi = trsidrji) 

= AijXjdt + eaij(sj)dt + 

d~i^ij^mj ^ij'^mj)dWrfi (20a) 

dixips) = {dxi)ps + Xidps + dxidps 
= A^jXjpsdt + eaij{sj)psdt + 

A{^ijC'mj ^ij'^mj)i‘^CmkPkdt + psAW^n) 

—ixie[sj,pj\dt + 2xiCmjPjdWm (20b) 

d{xiXjps) = AikXkXjpsdt + XiXkPsAljdt + 

A{^if^Cmk ^ik'kklmk^i^ jl^ml ^jlkklml)psdt^ 

= AikXkXjpsdt + XiXkPsAljdt 

+ + 2N,,)psdt, (20c) 

where we used the Riccati equation (5b) in the last equa¬ 
tion; moreover, we used Pi = Pi — Xips = 0(e). Finally, 
we also dropped stochastic terms in the last equation 
since they would give rise to a stochastic contribution of 
second order in e in the effective stochastic master equa¬ 
tion. 

Eq. (20c) is a Lyapunov equation; generally the 
steady-state solution of a Lyapunov equation AX + 
XA^ + B = 0 can be written as 

poo 

x = (21) 

Jo 

A straightforward calculation shows that in this case this 
amounts to 

x,x,ps = l{ri-^)ps, (22) 

which can be plugged into Eq. (18), which thus gets the 
form 

1 

dpi = AijPjdt - -eF"- [sj,ps]dt -t -ea^j [sj, ps]+dt -h 

here we used p = ps ® Pt + 0{e) which leads to 
tvT{[ri,rjAp) — 2piXj = T^jPs- We formally solve this 
equation; a straightforward calculations leads to 

''P = 2^^d^^Jk[sk,Ps]dt- -eA:r.^ajk[sk,Ps]+dt- 

A^j (F^-^Cm/c ^jk'kkimk)psdAVm’ (^^) 

We can already see that the unconditional part of the re¬ 
duced equation will not depend on the conditional state, 
as expected; since pi enters Eq. (16) only in the stochas¬ 
tic term, the unconditional part of (24) gives the only 
contribution to the unconditional dynamics of the sys¬ 
tem density operator ps- 

We proceed similarly to obtain an equation of motion 
for Pi. Combining Eqs. (18), (20b) and keeping terms to 


first order in e, we have 
= dr]i - d{xips) 

— -^ijPjdt j^^CiYik ^ik'^7nk)C-mjP‘jdt -\~ 

1 Z 

where fly = Uij — FJ^pg. The quantities fly can be 
interpreted in a similar way as the pi in Eq. (17): fly 
measures the deviation of the exact state p from the ten¬ 
sor product state pr ® ps with respect to the second 
order moments of the transducer’s canonical variables . 
The equation of motion for fly can be derived in a sim¬ 
ilar way as for pi, and shows that this is a second order 
quantity fly = O(e^). Therefore, the stochastic term can 
be dropped in Eq. (25), and the solution is 

1 z 

pi = ~l^^Qij 'Xjk[sk — (Sfe),ps]+ -|- —cQij Fjj,[sfe,ps], 

(26) 

where Q = A — 2(F'’cm — cnnm)cl. Plugging the results 
(24), (26) into the equation of motion of the system den¬ 
sity operator, Eq. (16), we recover the effective equation 
(8). To show that the resulting equation is a valid Be- 
lavkin equation, we need to show that each measurement 
channel has a corresponding decay term. This issue is 
addressed in Appendix B. 


C. Oscillating system operators 

When moving to the rotating frame with respect to 
the system Hamiltonian, the interaction stays time in¬ 
dependent only for special cases. Generally, the system 
operators will become time dependent. To go beyond the 
model presented in Sec. IIB, we now consider the sim¬ 
plest case of time dependent operators—those oscillating 
at frequency ±a;. We thus write the interaction Hamilto¬ 
nian as Hint = where s{t) = s+e“‘-|-s_e““‘, and 

the operators s± are time independent. Although this is 
not a completely general form of system-transducer cou¬ 
pling, together with the time-independent case, it can 
cover a large range of scenarios, including arbitrary qubit 
dynamics. 

Since the system operators now oscillate at frequency 
w, the essential part of the signal will no longer be trans¬ 
mitted in the carrier frequency of the transducer but by 
the sidebands instead. To recover this signal, we perform 
the measurements with local oscillators that are detuned 
from the standard reference frame. Denoting the fre¬ 
quency of the standard reference frame (corresponding, 
e.g., to the frequency of the laser light used for the read¬ 
out) as Wq and the frequency of the local oscillator as 
we can follow the approach of Ref. [33]. This adjust¬ 
ment results in time-dependent measurement operators 
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Am = {Cm + *Wm)^?'(i), where we have 






. aJe'A"** — a,;e 


— iAmt 


Qi = 


V2 


Pi = *- 


V2 


(27) 


with Am = Wo — Wm- Alternatively, we can also 
rewrite the measurement operators so that the time 
dependence enters through the coefficients, Xm{t) = 
{cm{t) + which will prove useful when adi- 

abatically eliminating the transducer dynamics. Overall, 
the stochastic master equation thus takes the form 


dp = -ie[s'^{t)r,p]dt + Crpdt + '^TL[Xmit)]pdWm, 

m 

(28) 

where we explicitly write the time dependence of the in¬ 
teraction Hamiltonian and the measurement operators. 

Before proceeding with the elimination procedure, 
some attention has to be paid to the conditional steady 
state of the Gaussian system. Since the measurement 
terms are now time-dependent, the Riccati equation (5b) 
for this system is ill-defined. To circumvent this problem, 
we perform a rotating wave approximation in the mea¬ 
surement terms by introducing the coarse-grained Wiener 
increments 


Since we made no assumptions about time dependence 
of the system operators in deriving equations of motion 
for PS, Pi, XiXjp, and pi, Eqs. (16), (18), (20c), (25), 
these equations are valid also in the present case. It is 
only their formal solution, where the time dependence of 
the system and measurement operators starts to play a 
role. The solution is, nevertheless, analogous to the time 
independent case, only with additional oscillation terms. 
Solving the equations of motion for XiXjps, Pi, 
Pi formally and performing the rotating wave approxi¬ 
mation, keeping only stationary terms, a straightforward 
calculation recovers Eq. (13). To bring the deterministic 
part of this equation to Lindblad form, we can proceed 
similar to the time-independent case. Since now the sys¬ 
tem operators are non-Hermitian, we first need to 
express them using some Hermitian basis (in the case of 
qubits, for instance, that would be the set of the Pauli 
operators and the identity). We can then recover the 
Hamiltonian part and the dissipative part, the diagonal- 
ization of which reveals the individual decay channels. It 
then remains to show that the measurement channels are 
included in the decay, for which we refer to Appendix B. 


III. EXAMPLES 


dW^ = j V2cosiAmt)dW, (29a) 

dW^ = J V2sm{A^t)dW. (29b) 

Eor integration intervals long on the time scale of 
but short on all other time scales, this produces two in¬ 
dependent Wiener increments, dW^dW^ = SmnSabdt, 
a,b = {c,s}, effectively turning every measurement into 
two, 

mXmrnpdw^ ^ -^n[xi,]pdw^ + -^n[xt^]pdw^, 

(30) 

where X!^ = (c((j -I- r and 

Cm{t) = cos(Amt) -h sin(Amt), (31a) 

mrait) = rrf^ cos(Am,t) -|- sin(Amf). (31b) 

These measurement operators are time-independent and 
thus give rise to a valid Riccati equation 

rc = AT^ + T'^A^ + 2N- (32) 

m aG{c,s} 

We treat the matrix Q [Eq. (9b)] which now also becomes 
time-dependent, in a similar manner; it becomes 

Q = (33) 

m a 

With these adjustments, we are now ready to adiabati- 
cally eliminate the transducer dynamics, and obtain an 
effective equation for the system. 


In this section, we illustrate the use of the adiabatic 
elimination method presented in Sec. H on a few simple 
examples. The model scenarios are taken from circuit 
QED where thermal noise—typically not accounted for 
by other adiabatic elimination methods—can be present 
even in cryogenically cooled systems. We show that our 
adiabatic elimination method, which we henceforth refer 
to as Gaussian adiabatic elimination, can provide signif¬ 
icantly increased accuracy for thermal noise at the level 
of few quanta. 

The examples we consider are illustrated in Eig. 2. 
In Sec. HI A, we consider dispersive readout of a qubit 
from a cavity that is coupled to a thermal reservoir, see 
Eig. 2(a). We compare Gaussian adiabatic elimination 
with results obtained by density operator expansion and 
show that significant qualitative and quantitative im¬ 
provements can be achieved with the former method. 
We extend this system in Sec. HIB where we study the 
effect thermal noise has on generating two-qubit entan¬ 
glement by measurement, following the approach of Ref. 
[19] [Fig. 2(b)]. Next, we illustrate the use of Gaussian 
adiabatic elimination with time-dependent interaction in 
Sec. HIC where we consider a single qubit coupled to 
a cavity field via Jaynes-Cummings Hamiltonian. Fi¬ 
nally, in Sec. HID, we consider the system shown in Fig. 
2 (c) —a transducer consisting of two coupled oscillators, 
one of which is coupled to a thermal bath. This setup 
differs from all other scenarios considered here by having 
a different unconditional and conditional steady state of 
the transducer, and we show how the Gaussian adiabatic 
elimination fares in this case. All numerical calculations 
in this section are done using QuTiP [34, 35]; for more 
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Figure 2. Schematic illustrations of the setups we consider 
to illustrate the Gaussian adiabatic elimination method. In 
Secs. Ill A, me, we analyze dynamics of a qubit coupled to 
a thermal cavity via dispersive and Jaynes-Cummings inter¬ 
action, respectively, shown in (a), (b) Setup for entanglement 
generation by measurement as discussed in Sec. IIIB. Here, 
two qubits interact dispersively with the same cavity but not 
with each other. Using the measurement record, it is then 
possible to postselect an entangled state of the two qubits. In 
Sec. HID, we study a qubit coupled to a two-oscillator trans¬ 
ducer, where the first oscillator couples to a thermal bath and 
the second oscillator is used for readout of the qubit state, as 
shown in (c). 


details on the particular implementations of these exam¬ 
ples, the reader ought to refer to Ref. [36]. 


A. Single qubit QND readout 

We consider the system shown in Fig. 2(a), where a 
qubit couples in a QND interaction to a cavity mode 
whose output hie is subject to continuous homodyne de¬ 
tection. In such a system the cavity itself serves just 
as a transducer, and can be adiabatically eliminated if 
the cavity decay rate is sufficiently large (faster than the 
QND coupling). In this case, adiabatic elimination is 
usually based on expanding the density operator in the 
Fock basis of the cavity around its vacuum state assum¬ 
ing no thermal excitations in the cavity [17, 19], 

P = Poo|0)(0l-|-pio|l)(0l-|-poi|0)(ll+piill)(ll + ..., (34) 

where the elements ptj are operators acting on the Hilbert 
space of the system, and are of the order i+j in the small 
coupling parameter e. Expanding up to second order, 
the reduced state of the qubit is given by ps = try (p) = 
Poo + Piij and the elements poO; Pii depend on p^ with 
i + j < 2. However, such an approach is limited to zero 
temperature where the cavity is essentially in vacuum. 
Our method allows to drop this assumption, and take 
thermal noise in the cavity into account in a systematic 
fashion. 

Before we illustrate our method on the basis of this 
example we note that more refined versions of adiabatic 
eliminations exist which employ a polaron-like transfor¬ 
mation [20, 23], and cover regimes of strong interactions 
between cavity and qubit. We believe that similar ap¬ 
proach, i.e., using different conditional steady states for 


different states of the system, is possible also with Gaus¬ 
sian adiabatic elimination; such a generalization is, how¬ 
ever, beyond the scope of the present paper. 

We start from the standard dispersive interaction with 
the qubit-cavity Hamiltonian of the form 

H = —Uz + Aa^a J- ga^aaz + (e*a + ea^). (35) 


We move to the interaction picture with respect to the 
free qubit evolution, canceling the first term. The second 
term gives the free cavity dynamics; we choose to drive 
the cavity mode at the center frequency, A = 0, maxi¬ 
mizing the measurement efficiency. The third term gives 
the standard dispersive interaction, and the last term de¬ 
scribes the cavity drive. 

To obtain an interaction that is linear in the cavity 
quadrature operators, we linearize the Hamiltonian by 
moving to the displaced frame, p —>■ {a)pD{a) with 

D{a) = exp(Q:a^ — a*a) being the displacement oper¬ 
ator, and a = —2ieln^ where k is the cavity decay 
rate. (The linearization also makes it possible to elim¬ 
inate the cavity field using density operator expansion 
approach.) This procedure brings the interaction Hamil¬ 
tonian to the form g(a*a+aa^)az+ga^aaz- If the driving 
field is strong enough, we can drop the second term, get¬ 
ting the interaction Hamiltonian iJint = where 

+ a'l’e"®'^) and y = y/2g\a\. The phase (j) is 
set by the field e driving the cavity. 

Since the cavity field couples to a thermal bath, the 
measurement term takes the form “i" ~ 

na^]p [1]. The full dynamics of the qubit-cavity system 
is thus described by the equation 


dp = —ixicTzVcf,, p]dt + K{{n + l)'D[a] J- nT>[ai]}pdt -I- 
-'H[{n + \)a — fia^pdW. (36) 


2ti - j -1 


Here, we assume that the cavity leaks only through its 
output port at rate k and the homodyne detector has 
unit efficiency; our numerical simulations indicate that 
additional decay has little effect on the accuracy of the 
adiabatic elimination methods. The measurement signal 
has the form dl = y/2Kl (2h -I- l){q)dt + dW. 

Following the recipe from Sec. HA, we have for the 
transducer Hamiltonian Ht = 0, jump operators ji = 
\lK{n + 1) a, j2 = and measurement operator 

A = Kj{2n + l)((h + l)a — ha^), or 


R = 0, 




niji + 1 ) 



(37a) 

(37b) 


^2 = \ ^ 


Kn I I 

^ V 

K 

2{2n + l) 

K{2n 1 ) 


(37c) 

(37d) 

(37e) 




















It then follows that A = —§1, JV = (h + i)l (1 be¬ 
ing the identity matrix), and both the unconditional and 
conditional steady state is the thermal state = 

(2fi -I- 1)1. Furthermore, from the interaction Hamilto¬ 
nian, we can read off s = X'^z(cos(f>, — sint/))^. Plugging 
everything into Eq. (8), a straightforward calculation 
reveals the effective equation 


2y^ 

dps = — (2n -b l)V[a^]psdt + 

K 


(38a) 


2x2 


K,(2n -b 1) 


% [—i{2ncos(j) + e *'^)crz] psdVF, 


d/ = —-1 


8x" 


K{2n -b 1) 


sm(l){az)dt + dW. 


(38b) 


Obviously, the optimal phase for an efficient readout of 
the qubit state is (j) = 7’'/2, which is not surprising—this 
phase choice corresponds to an interaction of the from 
Hint = X^zP accompanied by a g measurement. 

In contrast, using the density matrix expansion 
method, the qubit equation of motion takes the form 

2y2 

dps = -^{2,n + l)F>[cr^]psdf -b 


d/ = - 


+ \l^n[e-^^^+^/^^az\psdW, (39a) 

V K 

sin4i{(Tz)dt + dW. (39b) 


Sx' . 


Here we took into account the effect of thermal noise in 
the deterministic part (first line), which can be easily 
done using, e.g., projection operator method [37]. The 
only difference between Eqs. (38) and (39) is thus in the 
measurement term. Qualitatively speaking, the density 
operator expansion approach overestimates the strength 
of the measurement by a factor of ~ 1 /\/%. This means 
that for a zero temperature bath, both methods give the 
same results. In the presence of thermal excitations, how¬ 
ever, this difference quickly starts to play a role. 

To quantify the difference between the full model given 
by Eq. (36) and the effective qubit equation (38) or (39), 
we calculate the trace distance between the correspond¬ 
ing qubit states (we use pi to denote state obtained from 
the exact dynamics and p2 for the approximation meth¬ 
ods), D{pi,p 2 ) = ^trjpi — P2I with jXj = VX^X. Since 
the density matrices describe the state of a single qubit, 
the trace distance can be expressed using the expectation 
values of the Pauli matrices (aj) = tr(pj crQ as 


D{pi,P2) = 2 



(40) 


To obtain an average trace distance between the full 
model and the reduced dynamics, we generate a large 
number of quantum trajectories. We are thus able to 
study how the average trace distance changes in time; in 
addition, upon time averaging, we obtain a single figure 




20 30 

Time (k'') 



Figure 3. Illustration of determination of the average trace 
distance. Starting from a single quantum trajectory [expec¬ 
tation values of the Pauli matrices shown in plots (a)-(c)], we 
calculate the trace distance between the full model, Eq. (36), 
(dotted blue line) and the result obtained by Gaussian adia¬ 
batic elimination, Eq. (38), (full green line) or using the den¬ 
sity operator expansion, Eq. (39), (dot-dashed red line). The 
resulting trace distances are shown in (d). We further average 
using 500 quantum trajectories (e) to obtain an average trace 
distance. Using time averaging on this result, we further ob¬ 
tain a single figure of merit that determines quality of the two 
approaches. For the results shown here with n = 2, x ~ O.lre, 
0 = 7r/2 and initial qubit state Itpo) = (|0) -b \ V})/t/2, we have 
the average trace distance D « 0.05 for Gaussian adiabatic 
elimination and D « 0.22 for density operator expansion (cf. 
Fig. 4). 


of merit quantifying the discrepancy between the full and 
reduced dynamics; the averaging process is illustrated in 
Fig. 3. 

The results of the numerical investigations are shown 
in Fig. 4. In (a), we plot the average trace distance as 
a function of the interaction phase <(> for Gaussian adi¬ 
abatic elimination (green squares) and density operator 
expansion (black stars). We can see that both methods 
provide best results for </> = 7r/2, corresponding to an 
interaction of the form Hint = XP^z ■ This feature is par¬ 
ticularly beneficial since, as discussed before, this phase 
choice is optimal for non-demolition readout of the qubit 
state. 

In panel (b) of Fig. 4, we plot the average trace dis¬ 
tance versus thermal occupation number. While the av¬ 
erage trace distance with Gaussian adiabatic elimination 
(green squares for cj) = 7r/2, blue circles for tf) = 0) even¬ 
tually saturates (with the phase 4> = 0 this happens at 
h Ri 3, which is not shown in the plot), the error with 
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of view of solution accuracy. As there is no preferred 
phase for the qubit, the adiabatic elimination methods 
perform similarly for all these states. We also remark 
that generating quantum trajectories with the approx¬ 
imation methods is, due to smaller size of the Hilbert 
space, about four times faster than with the full model; 
in systems with larger thermal noise, this effect will be 
even larger. Moreover, as the qubit dynamics happens on 
a slower time scale than the evolution of the cavity field, 
it is possible to use larger time steps in the numerical 
solution, speeding the numerics up even more. 


B. Two-qubit entanglement by measurement 


Figure 4. Average trace distance for Gaussian adiabatic elim¬ 
ination and density operator expansion as compared to the 
full model. In (a), we plot the trace distance as a function of 
the measurement phase with green squares showing results for 
the Gaussian adiabatic elimination and black stars for den¬ 
sity operator expansion. The bottom panels show the trace 
distance versus thermal occupation number (b) and the over¬ 
all measurement time (c) for two choices of phase— (p = 7r/2 
(green squares for Gaussian adiabatic elimination, black stars 
for density operator expansion), and <j) = 0 (blue circles for 
Gaussian adiabatic elimination, red crosses for density oper¬ 
ator expansion). The parameters used for the simulations are 
X = O.lre, n — 2 [for (a) and (c)], measurement time Tm = 50 
[(a) and (b)], and initial qubit state \tpo) = (|0) -I- \l))/\/2. 
The Fock space of the cavity field for the full model is cut off 
at Amax = 20. 


density operator expansion (black stars for (/) = 7r/2, red 
crosses for </> = 0), as expected, grows with increasing 
temperature. Moreover, the Gaussian adiabatic elimi¬ 
nation performs a factor of about 2 better than density 
operator expansion already for half a thermal excitation 
present in the bath; with the phase choice (p = 7r/2, which 
corresponds to the optimal qubit readout, the difference 
between the two methods quickly grows. 

In Fig. 4(c), we investigate how the measurement time 
affects the accuracy of the two methods. Gaussian adi¬ 
abatic elimination remains unaffected by the length of 
the measurement {(p = 0) or even improves with time 
(</) = tt/2), whereas accuracy of the density operator ex¬ 
pansion method slowly deteriorates over time. This fea¬ 
ture can be seen already from the time-dependence of the 
trace distance [cf. Fig. 3 (e)], where the trace distance 
with Gaussian adiabatic elimination reaches a maximum 
shortly after the begin of the evolution {t « 5) and then 
settles at a smaller steady state value, while the trace dis¬ 
tance with density operator expansion continues to grow 
throughout the evolution. 

Finally, we note that the choice of a single initial qubit 
state |'0o) = (|0) -f \\))/\/2 does not affect the complete¬ 
ness of our analysis. Since the evolution for eigenstates of 
the Uz operator is trivial, the dynamics starting from the 
eigenstates of ax,y is the most interesting from the point 


Extending the system presented in the previous sec¬ 
tion, we now consider two qubits dispersively coupled to 
a common cavity field, Hint = 'Yhj where acts 

on the jth qubit. Such a system is of particular interest as 
the joint measurement of the two qubits can generate en¬ 
tanglement between them, as discussed in [19], recently 
realized experimentally [10] in a similar scenario. In¬ 
deed, a straightforward generalization of Eq. (38) (with 
(p = —7r/2) gives the effective dynamics 


dps = -(2n -I- 1 ) 1 ) 

K 


E 


k(2u + 1) 


n 






Psdt + 

ai psdFF, (41a) 


= + (41b) 


We thus get an effective measurement of the number of 
excitations of the two qubits. If we now prepare the 
qubits in the state Ji/'o) = ^(|0) + |1))G(|0) + 11)), engineer 
the interactions so that Xi = X 2 = Xi ^nd postselect only 
those trajectories with measurement o']. + a1 ~ 0, the 
two-qubit state takes the form j4'+) = (jOl) -b \lD))/\/2 
since there is one excitation in the system but we have 
no information on which of the two qubits is excited. 
Moreover, this state is also a dark state of the Lindblad 
term Via] -b o‘l]pg so it is a conditional steady state of 
the stochastic master equation (41). 

We note that this approach requires post-selection and 
thus generates entanglement only probabilistically. Using 
the dispersive interaction in the form Hi„t = goP<i{o] -b 
tr^), one can achieve also a true parity measurement cr^cr^ 
where the initial state as above collapses either onto ['1'+) 
or = (jOO) -b \ll))/\/2, generating entanglement be¬ 
tween the two qubits deterministically [9, 23, 38]. Al¬ 
though we believe it possible to generalize our Gaus¬ 
sian adiabatic elimination to include also coupling that 
is quadratic in the quadrature operators, we leave such 
analysis for future work and focus here on the probabilis¬ 
tic protocol only. 
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Figure 5. Example histograms of the integrated current at 
the beginning of the readout [t = (a)] and at a later 

time \t — 100k~^, (b)]. In panel (c), we plot the probability 
of the two qubits to be in the state | 00 ) (green dot-dashed 
line), |11) (red dotted line), and |'I'+) (blue full line). The 
simulations were run with the parameters x = 0.1k, h = 
0 , and a thousand trajectories were used for generating the 
histograms. 

In more detail, the entanglement is generated using 
the following protocol: First, the qubits interact with the 
cavity mode (we assume Xi = X 2 = x) and the output 
field is measured which is described by Eq. (41). After a 
time Tm, we have the integrated current 

J[T^) = / ™ dt/(t); (42) 

Jo 

if the integrated current is close to zero, the expecta¬ 
tion value {(t], -I- a1) = 0 and the qubits are in the state 
Idi-i.). The whole procedure is illustrated in Fig. 5. At 
an early time in the evolution [panel (a)], the distribu¬ 
tion of the integrated current is Gaussian but at a later 
time [panel (b)] three distinct peaks form with the cen¬ 
ter one corresponding to the qubits in the state |4l+). 
Quantitatively, the postselection is performed by using a 
threshold v and keeping the state iff |J| < z/. A small 
threshold thus results in a pure entangled state, albeit 
with a small success probability; increasing the thresh¬ 
old value, in turn, results in a mixed state with reduced 
amount of entanglement. 

We plot the results of the numerical simulations in 
Fig. 6 . We analyze the logarithmic negativity [39] of 
the resulting postselected state (full blue line) and the 
corresponding success probability (dashed green line) for 
cavity coupled to a vacuum bath [panel (a)] and a bath 
with h = 2 (c). Generally, in the presence of thermal 
photons, longer measurement times are needed to reach 
a maximally entangled state—we use the measurement 
times Tm = 100 k“^ for zero temperature bath in panel 


Figure 6. Logarithmic negativity (full blue line) and suc¬ 
cess probability (dashed green line) versus the postselection 
threshold v for h = 0 (a) and h = 2 (c). The measurement 
time is Tm = 100k“^ in (a) and Tm = 250k“^ in (c); more¬ 
over, in the insets, we plot the logarithmic negativity and 
success probability for Tm = 75k~^ (left) and Tm = 150k“^ 
(right) for both (a), (c). In addition, in panels (b), (d), his¬ 
tograms of the integrated currents corresponding to the re¬ 
sults in (a), (c) are shown. We use the coupling y = 0.1k and 
average over 1000 quantum trajectories. 

(a), and Tm = 250k“^ for the results plotted in panel 
(c). This effect is due to spreading of the peaks in the in¬ 
tegrated current with growing environment temperature, 
cf. Fig. 6 (b), (d). There, one can see that the local min¬ 
ima between peaks are slightly less pronounced for n = 2 
even with a measurement that is longer by a factor of 
2.5. 

Our observations are further accentuated in the insets 
of Fig. 6 (a), (c), where we plot the logarithmic negativ¬ 
ity and success probability for Tm = 75 k~^ (left inset) 
and Tm = 150k“^ (right inset). With thermal photons 
present, the logarithmic negativity does not reach unity 
in the limit zz —>■ 0 for the shorter time while with vacuum 
bath, a plateau of unit entanglement starts to form. For 
longer time, we reach a large plateau of success proba¬ 
bility of 0.5 with zero temperature, making it possible to 
generate the |'I'+) Bell state in half the cases; a similar 
plateau with the thermal bath starts to form only around 
Tm = 250 k-\ 


C. Single qubit Jaynes-Cummings readout 

To illustrate the adiabatic elimination with oscillating 
system operators, we consider a simple example of a sin¬ 
gle qubit coupled to a cavity mode via Jaynes-Gummings 
Hamiltonian, H = -|- Aa^a -|- ^(acr+e®'^ -I- aV_e“®'^). 

Moving to the rotating frame of the qubit, this gives rise 
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to the stochastic master equation 

dp = + Aa^a, p]dt + 

+K{n + l) 2 ?[a]pdt + Kh 2 ?[a^]pdt + 

+ Y/ ^_^ ^'H[(h + l)ae“*'*‘ — rla^e*'**]pdhh. (43) 

To obtain a full model without oscillating measurement 
operators, we now move to the rotating frame of the cav¬ 
ity. Choosing A = w = —6, Eq. (43) simplifies to 

dp = —tp[atT+e*‘^ -I- p]dt -I- 

+K.{n + l)2?[a]pdf -I- KfiT>[a^]pdt -I- 

-I--I-l)a — ha’^JpdlT, (44a) 

dl=^^^{q)dt + dw. (44b) 

Eq. (44) will be used in numerical calculations for com¬ 
parison with the adiabatic elimination methods; it is Eq. 
(43), however, that will be used as a starting point for 
the elimination of the cavity dynamics. This choice en¬ 
ables us, in principle, to go beyond the scenario with 
A = w = — (5 in the adiabatic approximation—using the 
Gaussian adiabatic elimination method, it is possible, for 
instance, to describe dynamics with measurement per¬ 
formed at the other sideband, S = oj. 

The transducer dynamics is given by the Hamilto¬ 
nian Ht = Aa^a, jump operators ji = y/K(n -\- 1 ) a, 
72 = y/1^ , and measurement operator A = 

y/K/{2n + l)((h-|-l)ae“*'^*—ha’I'e*'^*), so we have R = Al, 
same as for the dispersive readout in Eqs. (37b), 
(37c), and the measurement 


A = 


2{2fi+l) 


{cos((5t) — i(2n-I-1) sin((5t)}g (45) 


-h 


2 ( 2 rl-H 1 ) 


{sin((5<) -|- i{2n + 1) cos{St)}p, 


c = 


c = 


c+ = 


2{2n+l) \ 0 


2{2n+l) \ 1 



1 K{2n + 1 ) / 

), m 

' 2 ( 
1 K{2n + 1 ) / 

),rn^ = \ 

2 ( 


-1 


8(2h-h 1) 


1 

—i 


m~^ = 


K{2n + 1) f i 
1 


and c = (c+)*, m = (m+)*. We thus have A = —f 1 + 
Act, N = (n + i)l, and the cavity steady state (both 
unconditional and conditional) is the thermal state T'^ = 
T'^ = (2h -1-1)1. Together with the system operators 


^ A'**"- 


1 

—i 


(46a) 

(46b) 



Figure 7. (a) Average trace distance for the Gaussian adi¬ 
abatic elimination (green squares) and density operator ex¬ 
pansion (black stars) as a function of the interaction phase (p. 
In the bottom panels, we plot the average trace distance ver¬ 
sus thermal occupation (b) and measurement time (c) for the 
choice of phase 0 = 0 (blue circles showing Gaussian adiabatic 
elimination and red crossess for density operator expansion) 
and <j) = TV j2 (Gaussian adiabatic elimination shown in green 
squares, density operator expansion in black stars). The pa¬ 
rameters used to run the simulations are g = 0.1k, h = 2 [for 
panels (a), (c)], measurement time Tm = 50 [for (a), (b)], and 
initial qubit state |V^o) = (|0) -I- |l))/%/2. The cavity field for 
the full model has been cut oft at the Fock number Amax = 20. 


and the choice of frequencies A = ui = —S, this gives the 
stochastic master equation 

4(7^ 

dps = —{(n -b l)'D[a-] + n'D[a+]}psdt + 

Hi 

+ , n[{n + l)a_e-*(^+"/2) - 

y/Hi{2n + 1 ) 

—(0+^/2)]^dW, (47a) 

dl = {ay cos 0 — (Ta; sin 0 )dt + dH(.47b) 

y/ K{2n + 1) 


Using the density operator expansion method, together 
with a correction for thermal noise in the Lindblad terms, 
the qubit dynamics is described by the equation 


4q^ 

dps =- {{n -b l)T>[tT_] -b n'D[a+]}psdt + 

K 

-b^'H[a_e-*('^+’^/2)]psdW, 

y/K 

dl = {ay cos (j) — ax sin 0)dt -b dW. 

Y Hi 


(48a) 

(48b) 


Both adiabatic elimination methods, Eq. (47), Eq. (48), 
give identical results for zero-temperature cavity bath. 

The average trace distance for the Gaussian adiabatic 
elimination and the density operator expansion is ana¬ 
lyzed in complete analogy with the dispersive readout in 
Fig. 7. The error is minimized for phase 0 = 0 (a), 
which corresponds to a ay measurement, while for a ax 
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measurement (phase (j) = 7r/2), the average trace dis¬ 
tance reaches its maximum. We note, however, that the 
Jaynes-Cummings readout is much less phase-sensitive 
than the dispersive readot scheme. Performance of the 
Gaussian adiabatic elimination does not depend on the 
thermal occupation [panel (b)] while the average trace 
distance with the density operator expansion increases 
as expected. Finally, for long measurement times [see 
panel (c)], the average trace distance for both methods 
gradually decreases as the measurement approaches a 
projective readout and the qubit approaches one of its 
conditional steady states | 0 ), | 1 ). 

D. Two-oscillator transducer 

All examples considered so far had one special property 
in common—the unconditional and conditional steady 
states of the transducer were equal. To show how Gaus¬ 
sian adiabatic elimination can be applied to systems 
where this is not the case, we now consider the following 
example, see Fig. 2(c): A qubit, our system of interest, 
couples to a harmonic oscillator by means of a quan¬ 
tum non-demolition interaction similar to the example 
in Sec. Ill A. This oscillator decays into a thermal bath 
and, at the same time, couples to another oscillator of 
much higher frequency so its reservoir is effectively in 
the ground state. Finally, we measure the output of the 
second oscillator. 

The density operator of the overall system has the form 

dp = —i[x(^zr4, + uJia^a + L02b^b -\- gqiq2,p\dt -|- 
-|- 7 (h -I- \)'D[a\p<lt + xnD[a^]p(lt -\- 
-\-KD\b\pdt + \fK'H^e''^]pdW, (49a) 

d/ = + 6^e-*‘^)dt -f dW. (49b) 

Here, a describes the first (i.e., thermal) oscillator while 
b is used for the second readout oscillator, and r^j, denotes 
a general quadrature operator of the thermal oscillator. 
Such a system can be realized by coupling superconduct¬ 
ing qubit to a mechanical oscillator [40, 41] and read¬ 
ing out the signal in the mechanical oscillator optically. 
Since the oscillator coupling has the form of standard 
linearized optomechanical interaction, iJosc = P 9 i<? 2 , the 
qubit readout is optimized by driving the readout oscil¬ 
lator on the red sideband, wi = W 2 = w. The readout 
efficiency can further be maximized by letting the qubit 
couple to the phase quadrature of the thermal oscillator 
and measuring the phase quadrature of the readout os¬ 
cillator. The stochastic master equation then takes the 
form 

dp = —ilxcTzPi + u}{a^a + b^b) + gqiq2,p\dt + 

-l- 7 (h -I- l)'D[a\pdt + xnDla^jpdt + 

+K'D[b]pdt + v^'H[i 6 ]pdIF. (50) 

As the transducer dynamics is more complex in this 
case, we perform the whole adiabatic elimination numer- 



Thermal occupation Time(K ') 

Figure 8. (a) Average trace distance as a function of thermal 
occupation for three regimes: weak coupling {g = 0.2k, green 
squares), intermediate coupling {g = 0.5k, blue circles), and 
strong couping {g = k, black stars). In (b), we show the 
average trace distance versus time for n = 2, g = k. Other 
parameters used in the simulations are y = 0.2k, lj = 5k, 
7 = 0.1k, initial qubit state Itpo) = (jO) -|- \l))/y/2, and we 
averaged over 100 quantum trajectories. 

ically, see Ref. [36]. The results of the numerical sim¬ 
ulations are shown in Fig. 8 . In panel (a), we inves¬ 
tigate how the bath temperature for the thermal oscil¬ 
lator affects the average trace distance in three distinct 
regimes—for weak {g = 0.2k, green squares), intermedi¬ 
ate [g = 0.5k, blue circles), and strong [g = k, black 
stars) coupling. As the strength of the coupling between 
the two oscillators grows, the trace distance becomes less 
temperature sensitive. In Fig. 8 (b), we plot an example 
average trace distance as a function of time. This plot il¬ 
lustrates that the time dependence has features similar to 
simpler transducers considered in the previous sections— 
after a short initial transient time, the trace distance sat¬ 
urates and stays constant for the rest of the evolution. 

Finally, we also point out the numerical requirements 
for the full model and the adiabatic elimination. With 
only two thermal excitations in the heat bath, the full 
model needs 700 times longer time to be simulated, com¬ 
pared to the adiabatic elimination; this difference can be 
increased by using larger time steps for the approximate 
dynamics since the qubit evolution happens at longer 
time scales. The main limitation in our numerical anal¬ 
ysis, however, are the memory requirements. With two 
thermal excitations (and corresponding Fock space cut¬ 
offs at 20 and 10 excitations for the thermal and read¬ 
out oscillator, respectively), the storing of the full density 
matrix for the whole time evolution requires several giga¬ 
bytes of working memory. Since the cutoff energy grows 
faster than linearly with increasing temperature, and the 
size of the density matrix grows quadratically with the 
cutoff, we were not able to perform reliable numerical 
simulations for larger bath temperatures. We still be¬ 
lieve, nevertheless, that Gaussian adiabatic elimination 
can be used for systems with tens or hundreds thermal 
excitations present. 

IV. CONCLUSIONS 

In summary, we presented a new method of adiabatic 
elimination of fast degrees of freedom from stochastic 
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quantum dynamics. Assuming the transducer (i.e., the 
system we wish to eliminate) is Gaussian, we can fully 
describe its evolution using first and second statistical 
moments of its canonical operators; moreover, the co- 
variance matrix of the conditional state obeys determin¬ 
istic Riccati equation. We are thus able to treat trans¬ 
ducers coupled to thermal bath or consisting of multiple 
modes. While eliminating several modes using the ap¬ 
proach based on density operator expansion or polaron 
transformation quickly becomes tedious, our method re¬ 
quires only basic linear-algebraic tools and can be easily 
solved numerically. 

Since the procedure we use relies on the fact that the 
system of interest itself has no free evolution, we did not 
present a completely general treatment—instead, we fo¬ 
cused on the most relevant situations only. In the first 
place, we assumed that moving to the rotating frame 
with respect to the free Hamiltonian of the system leaves 
its interaction with the transducer time independent, 
corresponding, in particular, also to a quantum non¬ 
demolition interaction. Secondly, we considered a sce¬ 
nario, where the interaction has terms oscillating at the 
frequency ±w. Adapting the method for other forms of 
coupling is straightforward. 

With these results, we have shown how our method 
can be used to simulate readout of qubits in cavity QED, 
as relevant e.g. to superconducting circuit QED. Com¬ 
pared with the method of expanding the density opera¬ 
tor around the vacuum state of the readout cavity, our 
method provides significantly better results already for 
a few thermal excitations present and is thus relevant 
to many experimental scenarios. We believe that fur¬ 
ther improvements can be achieved with ideas borrowed 
from adiabatic elimination using polaron transformation. 
There, for strong coupling between the system and the 
transducer, one has to consider different steady states of 
the transducer for individual states of the system and 
perform adiabatic elimination with respect to these con¬ 
ditional states. Using similar tools for our method, it 
should be possible to eliminate any Gaussian transducer 
with respect to several conditional steady states. In ad¬ 
dition, it is possible to generalize the method for system- 
transducer coupling that is quadratic in the canonical 
transducer operators—one could, e.g., use such a result 
to analyze a full dispersive readout of a superconducting 
qubit, i?i„t = ga^a'a. 

Another field that could benefit from our results is cav¬ 
ity optomechanics. Typical frequencies of mechanical os¬ 
cillations can correspond to thermal noise of a few hun¬ 
dred quanta even with cryogenic cooling. Such systems 
cannot be eliminated from stochastic master equations 
using present methods; a toy model in our last exam¬ 
ple shows how similar tasks can be achieved using our 
approach. In the future, it might be interesting to study 
how optomechanical systems used for conversion between 
microwave and optical fields [29, 30] could be used to en¬ 
tangle two superconducting qubits using measurement 
and feedback [42]. 
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Appendix A: Equations of motion for mean values 
and covariance matrix of a Gaussian system 

In this Appendix, we derive equations governing the 
evolution of the mean values and the covariance matrix 
of a Gaussian system. We start with the master equation 

p=-i[H,p\ + '^V[jn]p. (Al) 

n 

Since the system is Gaussian, the Hamiltonian is 
quadratic in the canonical operators and we can write 
H = Rr (linear Hamiltonian leads to a simple dis¬ 
placement which can be treated suitably moving the 
origin of the phase space); the jump operators are lin¬ 
ear, and we write Here r = 

{qi,Pi,..., qNjpjy)'^ is vector of the canonical operators 
whose commutation relations define the symplectic ma¬ 
trix 


1 

" 0 I\ 

/ 

0 1\ 

'lO'ij , (T — I 

-1 0 )®- 

..® 

-10 


The goal is to use the master equation to obtain equations 
of motion for the first and second statistical moments 
defined by 

X = (r) = tr(pr), T^ = ([ri,rj]+) - 2xiXj. (A3) 
For the i-th mean value, we have 

Xi = tv^pri) 

= p]r J + ^ ti{V[jn]pri} 

n 

= -iiv{p[ri,H]} + 

+ - ]^[ji3n,n]+)}. (A4) 

n 

The commutator in the first term can be rewritten as 

jk 

= Rjk{[ri,rj]rk + rjln, rfc]) 

jk 

= '^crijRjkrk, (A5) 

jk 
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where we used the fact that the Hamiltonian matrix is 
symmetric R = . For the Lindblad terms, we have 

ir{V[jn]pri} = '^UjCk^^{pirknrj - ^[rkrj,n]+)} 
jk 

jk 

~ ~2 iinj^nk ~ ^nj^nk)^k- (^ 6 ) 

jk 

Combining everything, we can write 

^ ' ^ijRjkXk 2 'y y ^iji^njink ^nj^nk^^ki (-^'^) 
j k nj k 

or, in the matrix form 

x = Ax, A = aR-^a^{^l^n-CC)- (AS) 

n 

For the covariance matrix, we need to evaluate 

tij = tr{p[ri,r-j]+} - 2{xiXj +XiXj). (A9) 

Similar to the previous case, we have for the coherent 
evolution 

[rirj,H] = Rki{ri[rj,rk]ri + rkin, ri]rj + 

kl 

+ [ri,rk]rjri + rfer^r^, n]), (AlO) 

which, combined with [rjri,H], gives 


where 

aki = lY.^Ck(ni + UCi) (AlSa) 

n 

Pkl = 2 ^i^nkih - Ck^nl)- (A15b) 

n 

Plugging everything into Eq. (A9) and using XiXj = 
-^ikXkXj, and writing the resulting expression in ma¬ 
trix form, we get the Lyapunov equation 

t = AT + TA^ + 27V, (A16a) 

(A16b) 

n 

When the dynamics is described by the conditional 
master equation 

dp = -i[H, p]dt + ^ 2?[j„]pdt -f ^ HiXmlpdWm, 

n m 

(A17) 

we also need to evaluate the contributions from the mea¬ 
surement terms 'H[Am]p- We start by splitting the mea¬ 
surement operator in its Hermitian and anti-Hermitian 
part. Am = (cm + irUmVr = X)fc(cmfc + immk)rk, so we 
can write 

'H[Xm]p = [c^{r - x), p]+ + i[m^r, p]. (A18) 

For the mean values, this gives the contribution 


[[ri,rj] + ,H] = i'^{aikRki[rj,ri]+ - [ri,ri]+Rikcrkj). 
kl 

(All) 

For the decay terms, we have 

£[jn]rirj =jlnrjjn - ^[jijuAirjU 

= Ck^^iiirk, nrjjn + rk [nr ^, n]) 

kl 

= C^k^micyjirkn + crurkTj - 

kl 

-o-jknn - aikVjri), (A12) 

where we used 

[riTj , Tfc] = ri [rj , r^] + [r^, 

= icTjkri -b iaikTj. (A13) 

Combined with £[jn]rjri and summed over n, this ex¬ 
pression gives 

'^^[jn][ri,rj]+ = '^{l3ki{(Jjk[ri,ri]+ + aik[rj,ri] + ) - 
n kl 

-2akicrjkcru}, (A14) 


dxi = tr('H[Am]p)dVFm 

= ^tr{cmfep([ri,rfc]+ - 2xkrip) + 
k 

Mmmkp[ri,rk]}dWm 

— ^ ikC-mk 
k 

The mean values thus obey the equation 

da; = Axdt-b ^^(Fcm — o'Wm)dlFm. (A20) 

m 

For the covariance matrix, we need to evaluate 

dVij = tr{[ri,rj]+n[Xm]p)dWm - 

—2{{dxi)xj + Xidxj + da;idxj}, (A21) 

where we used the ltd rule d(Ay) = (dA)F -b AdF -b 
dAdF for the contribution of the mean values. In the 
following we concentrate on the stochastic contribution 
in the increments da;^ (second term on the r.h.s. of 
Eq. (A20)) as the deterministic contribution is trivial. 
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We start by considering 

tr([ri,rj]+dp) ='^Cmktr{p{{rk, [ri,rj]+]+ - 
k 

-2xk[ri,rj\+)} + 

+i ^ mmk^v{p[rirj + rfc]} 

k 

— 2 ^ ^ CjnkiJ^ik^j “t“ 
k 

-2'^mmk{x^(jjk + Xjdik). (A22) 

k 

In the first sum on the right hand side, we used the fact 

([^z; [^J;^fc]+]+) — ‘^(J^ijXk “1“ ^jkXi “1“ ^kiXj 2XiXjXk^- 

(A23) 

This can be seen by comparing the third derivative of 
the characteristic function from the definition x(r) = 
tT{D{r)p} (here D{r) is the displacement operator) with 
a general Gaussian characteristic function 

x(r) = exp{—iT-^crx — -r^tr^ror}. (A24) 

We further use 
{dxi)xj + Xidxj + dxidxj = 

— ^ ', Xi jkCmk ^jk'^mk')dWjri T 

k 

+ ^ XjiTikCmk - (Jtkmmk)dWm + (A25) 

k 

T ^ '^ {^ik^mk 'XikTTljrik^^jl^Tnl ^jl‘^Tnl)dt. 

k,l 

Combining everything, the stochastic contributions to 
the covariance matrix cancel out, and we are left with 
the term 

(TCm - (Tmm)(rcm “ arUm)^. 


The dynamics of the covariance matrix is thus given by 
the Riccati equation 

r = AT + r + 2 N - 2 y^{Tcm - <xmm){^Cm - crnim)'^ ■ 

m 

(A26) 


Appendix B: Positive-semidefiniteness of decay after 
subtracting the measurement channels 

We start by the observation that the overall decay E 
is positive. Using its definition, Eq. (11), together with 
the Lyapunov equation (7b), and the definitions (6), we 
can write 

P = A-^ + ^{Aa - a^A'^)^ 

= (Bl) 

i 

Next, we discuss the question whether the effective 
stochastic master equations in Eqs. (8) and (13) present 
valid Belavkin equation, that is, whether they generate 
completely positive maps. In order for this to be true 
each measurement channel has to have a corresponding 
decay process. Quantitatively, matrix P in Eq. (11) de¬ 
scribing all decay terms needs to be larger than the ma¬ 
trix AmAjyj characterizing all measurement channels. 
In other words 

p' = p-J2a^aI 

m 

has to be positive semidefinite. We did not prove this 
statement in the general case, but checked it for all of 
the cases treated in Sec. III. 
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